//============================================================================
// Name        : FabVoronoi.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <iterator>
#include <iomanip>
#include "Topol.h"
#include "FAtoms.h"
#include "HeaderTrj.h"
#include "ResidueCM.h"
#include "FVoronoi.h"
#include "FComms.h"
#include "Help.hpp"
#include "FstreamC.h"
#include "FstreamF.h"
#include "VoronoiPrint.h"
#include "trjInput.h"
#include "InputRadii.h"
#include "RadiiDBase.h"
#include "Myradius.h"

static int pGroup=-1; // From 0 to ngrps, ngrps no. of selected groups
static bool bPrintVols=FALSE;
static bool bPrintAreas=FALSE;
static bool bPrintDomain=FALSE;
static bool bDoCorrel=FALSE;
static bool bSpheres=false;


using namespace Topol_NS;
using namespace std;
int main(int argc, char ** argv)
{
	trj::trjInput MyInput(argc,argv);
    Parallel::comm=new Parallel::MyMPI(argc,argv);
    bool bHyd=true;
	bool bDel=true;
	bool bWriteFab=false;
	bool bIsrd=false;
	bool bTestVol=false;
	int narg=1;
	string filein,fileout="Voronoi.out",filepdb,filedomain,filedefdomain,fileirl,fileidb;

	string scratchdir;
	ifstream fpdb;
	ifstream *firl;
	ifstream * fidb;
	ifstream ftest,fdefdomain;
	Fstream * finx;
	ofstream fout;
	ofstream fdomain;
	string comm0(argv[0]);
	size_t mypos=comm0.find_last_of("/")+1;
	size_t length=comm0.size()-mypos;
	string command=comm0.substr(mypos,length);
	string errmsg;
	string Usage="Usage:\t"+ command + "\n";
	vector<string> use=MyInput.getUsage();
	for(unsigned int n=0;n<use.size();n++)
		Usage+=use[n];
	Usage+="\n\t Default values in square brackets []\n";
	try{
		if(argc == 1) throw Usage;
		if(int m=MyInput.bTest().size()) {
			errmsg=" Command(s) not found: ";
			for(unsigned int n=0;n<m;n++)
				errmsg+=MyInput.bTest()[n]+"  ";
			errmsg+="\n"+Usage;
			throw errmsg;
			}
		}
	catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}



	int nstart=0;
	int nend=-1;
	int nskip=1;
	// input read
	bool inputfile=false;
	try{
		if(!MyInput["-dcd"].empty()) {
			if(MyInput["-dcd"].size() != 2) throw string(" filename expected for " + MyInput["-dcd"][0] + " option ");
			filein=MyInput["-dcd"][1];
			ftest.open(filein.c_str(),ios::in);
			if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
			ftest.close();
			finx=new FstreamF(filein);
			inputfile=true;
		}
		if(!MyInput["-xtc"].empty()) {
			if(MyInput["-xtc"].size() != 2) throw string(" filename expected for " + MyInput["-xtc"][0] + " option ");
			filein=MyInput["-xtc"][1];
			ftest.open(filein.c_str(),ios::in);
			if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
			ftest.close();
			finx=new FstreamC(filein);
			inputfile=true;
		}
		if(!MyInput["-pdb"].empty()) {
			if(MyInput["-pdb"].size() != 2) throw string(" filename expected for " + MyInput["-pdb"][0] + " option ");
			filepdb=MyInput["-pdb"][1];
			fpdb.open(filepdb.c_str(),ios::in);
			if(!fpdb) throw string("\n Cannot open " + filepdb + "!!\n");
		}
		if(!MyInput["-o"].empty()) {
			if(MyInput["-o"].size() != 2) throw string(" filename expected for " + MyInput["-o"][0] + " option ");
			fileout=MyInput["-o"][1];
		}
		if(!MyInput["-b"].empty()) {
			if(MyInput["-b"].size() != 2) throw string(" Number of first frame needed for " + MyInput["-b"][0] + " option ");
			stringstream(MyInput["-b"][1])>> nstart;
		}
		if(!MyInput["-e"].empty()) {
			if(MyInput["-e"].size() != 2) throw string(" Number of last frame needed for " + MyInput["-e"][0] + " option ");
			stringstream(MyInput["-e"][1])>> nend;
		}
		if(!MyInput["-skip"].empty()) {
			if(MyInput["-skip"].size() != 2) throw string(" Number of skipped frames needed for " + MyInput["-skip"][0] + " option ");
			stringstream(MyInput["-skip"][1])>> nskip;
		}
		if(!MyInput["-pvol"].empty()) {
			if(MyInput["-pvol"].size() != 1) throw string(" No argument for " + MyInput["-pvol"][0] + " option ");
				bPrintVols=true;
		}
		if(!MyInput["-pgroup"].empty()) {
			if(MyInput["-pgroup"].size() != 2) throw string(" Group number needed for " + MyInput["-pgroup"][0] + " option ");
			stringstream(MyInput["-pgroup"][1])>> pGroup;
		}
		if(!MyInput["-parea"].empty()) {
			if(MyInput["-parea"].size() != 1) throw string(" No argument for " + MyInput["-parea"][0] + " option ");
				bPrintAreas=true;
		}
		if(!MyInput["-domain"].empty()) {
			if(MyInput["-domain"].size() != 1) throw string(" No argument for " + MyInput["-domain"][0] + " option ");
				bPrintDomain=true;
		}
		if(!MyInput["-correl"].empty()) {
			if(MyInput["-correl"].size() != 1) throw string(" No argument for " + MyInput["-correl"][0] + " option ");
				bDoCorrel=true;
		}
		if(!MyInput["-autodomain"].empty()) {
			if(MyInput["-autodomain"].size() == 2){
				filedomain=MyInput["-autodomain"][1];
				fdomain.open(filedomain.c_str(),ios::out);
				if(!fdomain) throw string("\n Cannot open " + filedomain + "!!\n");
			} else if(MyInput["-autodomain"].size() == 1) {
				filedomain="autodomain.pdb";
				fdomain.open(filedomain.c_str(),ios::out);
				if(!fdomain) throw string("\n Cannot open " + filedomain + "!!\n");
			}
			else
				throw string(" Too many arguments for " + MyInput["-autodomain"][0] + " Give one argument to pick your "
						"own output .pdb file ");
			bPrintDomain=true;
		    Topol::bAuto=true;
		}
		if(!MyInput["-defdomain"].empty()) {
			if(MyInput["-defdomain"].size() == 2){
				filedefdomain=MyInput["-defdomain"][1];
				fdefdomain.open(filedefdomain.c_str(),ios::in);
				if(!fdefdomain) throw string("\n Cannot open " + filedefdomain + "!!\n");
			}
			else throw string(" filename expected for " + MyInput["-defdomain"][0] + " option ");
		    Topol::bDefDomain=true;
		}
		if(!MyInput["-w"].empty()) {
			if(MyInput["-w"].size() != 1) throw string(" No argument for " + MyInput["-w"][0] + " option ");
				bSpheres=true;
		}
		if(!MyInput["-irl"].empty()) {
			if(MyInput["-irl"].size() != 2) throw string(" filename expected for " + MyInput["-irl"][0] + " option ");
			fileirl=MyInput["-irl"][1];
			firl=new ifstream;
			firl->open(fileirl.c_str(),ios::in);
			if(!(*firl)) throw string("\n Cannot open " + fileirl + "!!\n");
		}
		if(!MyInput["-idb"].empty()) {
			if(MyInput["-idb"].size() != 2) throw string(" filename expected for " + MyInput["-idn"][0] + " option ");
			fileidb=MyInput["-idb"][1];
			fidb=new ifstream;
			fidb->open(fileidb.c_str(),ios::in);
			if(!(*fidb)) throw string("\n Cannot open " + fileidb + "!!\n");
		}
		if(!MyInput["-nohyd"].empty()) {
			if(MyInput["-nohyd"].size() != 1) throw string(" No argument for " + MyInput["-nohyd"][0] + " option ");
				bHyd=false;
		}
		if(!MyInput["-hyd"].empty()) {
			if(MyInput["-hyd"].size() != 1) throw string(" No argument for " + MyInput["-hyd"][0] + " option ");
				bHyd=true;
		}
		if(!MyInput["-nodel"].empty()) {
			if(MyInput["-nodel"].size() != 1) throw string(" No argument for " + MyInput["-nodel"][0] + " option ");
				bDel=false;
		}
		if(!MyInput["-del"].empty()) {
			if(MyInput["-del"].size() != 1) throw string(" No argument for " + MyInput["-del"][0] + " option ");
				bDel=false;
		}
		if(!MyInput["-rd"].empty()) {
			if(MyInput["-rd"].size() != 1) throw string(" No argument for " + MyInput["-rd"][0] + " option ");
				bIsrd=true;
		}
		if(!MyInput["-nord"].empty()) {
			if(MyInput["-nord"].size() != 1) throw string(" No argument for " + MyInput["-nord"][0] + " option ");
				bIsrd=false;
		}
		if(!MyInput["-fab"].empty()) {
			if(MyInput["-fab"].size() != 1) throw string(" No argument for " + MyInput["-fab"][0] + " option ");
				bWriteFab=true;
		}
		if(!MyInput["-nofab"].empty()) {
			if(MyInput["-nofab"].size() != 1) throw string(" No argument for " + MyInput["-nofab"][0] + " option ");
			bWriteFab=false;
		}
		if(!MyInput["-test"].empty()) {
			if(MyInput["-test"].size() != 1) throw string(" No argument for " + MyInput["-test"][0] + " option ");
				bTestVol=true;
		}
		if(!MyInput["-help"].empty()) {
			if(MyInput["-help"].size() != 1) throw string(" No argument for " + MyInput["-help"][0] + " option ");
				throw Help;
		}
	}catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}
	VoronoiPrint::bPrintAreas=bPrintAreas;
	VoronoiPrint::bPrintVols=bPrintVols;
	VoronoiPrint::bPrintDomain=bPrintDomain;
    VoronoiPrint::pGroup=pGroup;
    Topol::fDomainPDB=new std::ofstream ;
    Topol::fDomainPDB=&fdomain;


	Fstream & fin=*finx;
	try{
		if(!inputfile) throw string("\nNo input file given will operate on the .pdb file\n");
	}
	catch(const string & s ){
		cout << s << endl;
		bIsrd=true; // Must do radical tessellation if no periodic conditions!
		Parallel::comm->Barrier();
	}
	if(bSpheres) bIsrd=true; // Must do radical tesselation if atomic spheres are use as initial
	// Voronoi cells

// read pdb file to construct topology
	vector<string> data;
	for(string str;getline(fpdb,str);){
		data.push_back(str);
	}
	if(!bHyd) Topol_NS::Myradius::setNOHYD();
	Topol MyTop;

// Define radii to use
	InputRadii IR(firl);
	RadiiDBase RDB(fidb);
	Myradius Radical(IR,RDB);
	cout << Radical("AOT","C2") << endl;

	// def domain
	if(Topol::bDefDomain) {
		string Dom0;
		int begin0,end0;
		vector<int> tmp(2);
		while(fdefdomain >> Dom0 >> tmp[0] >> tmp[1]) {
			MyTop.AddDomain(Dom0,tmp);
		}
	}

	MyTop(data,bIsrd);
	int natoms=MyTop.Size();
	int nresid=MyTop.ResSize();
	FAtoms atm;
	atm.setDim(natoms);

// Map atoms to residues to use in Voronoi. Eliminate hydrogens if needed
	vector<int> * resmap=new vector<int> [nresid];
    string * ResLabels=new string [nresid];
    string label;

    for(int i=0;i<nresid;i++){
    	ResLabels[i]=MyTop.ResInfo(i);
    }
    for(int i=0;i<natoms;i++){
    	int j=MyTop.ResIndex(i);

    	if(bHyd) resmap[j].push_back(i);
    	else {
    		label=MyTop.AtomInfo(i);
    		if(label.compare(0,1,"H",0,1) != 0) resmap[j].push_back(i);
    	}
    }

    ResidueCM MyCM;
    MyCM(natoms,nresid,resmap,ResLabels);

    unsigned int totaltypes=MyTop.getCidxSize();
    AtomIndex * cidx=new AtomIndex [totaltypes];
    for(unsigned int i=0;i<MyTop.getCidxSize();i++){
		const vector<int> & ty=MyTop.getCidx(i);
		cidx[i].set(ty.size(),&ty[0]);
	}
	MyCM.ResidueSelect(totaltypes,cidx);

// Setup Voronoi

	FVoronoi * vor=new FVoronoi(natoms,MyTop.getAtomName(),bHyd);

	vor->setTypes(MyTop.getmaxt(),cidx,MyTop.getTypeNames());
	vor->setResTypes(MyTop);

	Voronoi * RegVor=dynamic_cast<Voronoi *> (vor);
	ios::streampos len;
	HeaderTrj header;
// Read header of dcd file
	if(finx) {
		fin.seekg(0,"end");
		len=fin.tellg();
		fin.seekg(0,"beg");
		fin>>header;
		try{if(!header.check(MyTop.Size())) throw string("Number of atoms in the pdb and trajectory files does not match!");}
		catch(const string & s){cout << s<<endl;exit(-1);}
	}


	atm.setrd(MyTop.getrd());

// Read trajectory and compute Voronoi data

    Parallel::comm->Barrier();

    vector<MySubunit> MySub=MyTop.getMySub();
    Parallel::FComms Comms(fout,fileout,nstart,nend,header.getNFR(),nskip);
    if(!finx){
    	atm.pdb(data);
		vor->StartNP(atm);
		vor->getData(MySub);
		Comms.getStream() << *RegVor;
		if(bTestVol) vor->testVol();

    	Parallel::comm->Barrier();
    	Parallel::comm->~MyMPI();
    	cout << "\nProgram completed: Output data written to " + fileout << "\n\n";
        return 0;
    }
	nstart=Comms.getStart();
	nend=Comms.getEnd();
	int ntime=0;

	while(fin.tellg() < len){
		if(ntime < nstart) {fin+=atm;ntime++;continue;}
		if(nend != -1 && ntime > nend) break;
		fin >> atm;
		cout << fixed << setw(5) << "----> Step No. " << ntime <<"\n";
		if(bSpheres)
			vor->StartNP(static_cast<float> (ntime),atm);
		else
			vor->StartP(static_cast<float> (ntime),atm);
		vor->getData(MySub);
		if(bWriteFab)
			Comms.getStream() << *vor;
		else
			Comms.getStream() << *RegVor;

		if(bDoCorrel) vor->Stats(static_cast<double> (ntime));
		if(bTestVol) vor->testVol();
		ntime++;
	};
	if(bDoCorrel) vor->WriteStats(Parallel::comm);
	Comms.appendStreams();
	if(bDel) Comms.removeFiles();
	Parallel::comm->Barrier();
	Parallel::comm->~MyMPI();
	cout << "\nProgram completed: Output data written to " + fileout << "\n\n";
    return 0;
}   //end of main function
